{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Derivative Removal by Adiabatic Gate\n",
    "*Copyright (c) 2021 Institute for Quantum Computing, Baidu Inc. All Rights Reserved.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Outline\n",
    "This tutorial will demonstrate how to implement an X gate employing the DRAG (Derivative Reduction by Adiabatic Gate) technique using Quanlse. The outline of this tutorial is as follows:\n",
    "- Introduction\n",
    "- Preparation\n",
    "- Define the waveform for DRAG\n",
    "- Quanlse realization\n",
    "- Summary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "In superconducting circuits, one has to consider the leakage error due to the fact that superconducting circuits are not perfect two-level systems. For weakly anharmonic qubits, leakage into the third energy level takes the qubit out of the computational subspace. To overcome this issue, researchers proposed the DRAG procedure \\[1\\], which removes most of the leakage error by modifying the waveforms of the drive pulses."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preparation\n",
    "After you have successfully installed Quanlse, you could run the Quanlse program below following this tutorial. To run this particular tutorial, you would need to import the following packages from Quanlse and other commonly-used Python libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import numpy and scipy\n",
    "import numpy as np\n",
    "from scipy import integrate\n",
    "\n",
    "# Import the Hamiltonian module\n",
    "from Quanlse.Utils import Hamiltonian as qham\n",
    "\n",
    "# Import the function for calculating infidelity\n",
    "from Quanlse.Utils.Tools import unitaryInfidelity\n",
    "\n",
    "# Import related operators\n",
    "from Quanlse.Utils.Operator import driveX, driveY, number, duff\n",
    "\n",
    "# Import waveforms and functions used to process the waveforms' data\n",
    "from Quanlse.Utils.Waveforms import gaussian, dragY1\n",
    "from Quanlse.Utils.Waveforms import makeWaveData, waveFuncToSeq\n",
    "\n",
    "# Import simulator interface for Quanlse Cloud Service\n",
    "from Quanlse.remoteSimulator import remoteSimulatorRunHamiltonian as runHamiltonian\n",
    "\n",
    "# Import matplotlib for graphing purposes\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To use Quanlse Cloud Service, we need to acquire a token to get access to the cloud. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import Define class and set the token for cloud service\n",
    "# Please visit http://quantum-hub.baidu.com\n",
    "from Quanlse import Define\n",
    "Define.hubToken = \"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define the waveform for DRAG\n",
    "By performing a rotating wave approximation (RWA), the Hamiltonian in the rotating frame can be written as \\[2\\]: \n",
    "\n",
    "$$ \n",
    "\\hat H_R / \\hbar = \\delta_1 |1\\rangle \\langle 1|+\\delta_2 |2\\rangle \\langle 2|+\\alpha_q\\hat a^{\\dagger}\\hat a^{\\dagger}\\hat a \\hat a/2+\\frac{\\varepsilon_x(t)}{2}\n",
    "\\left[ \\hat{a}^\\dagger + \\hat{a} \\right]+\\frac{\\varepsilon_y(t)}{2}\n",
    "i \\left[\\hat{a}^\\dagger - \\hat{a}\\right]\n",
    ",\n",
    "$$\n",
    "\n",
    "where $\\omega_1$ and $\\omega_2$ are the qubits' frequencies; and $\\omega_d$ is the driving frequency. $\\alpha_q = \\omega_2 -2\\omega_1$ is the anharmonicity of the system. $\\delta_1 = \\omega_1-\\omega_d$ and $\\delta_2 = \\omega_2-\\omega_d$ are the detunings of the transitions with respect to the drive frequency. $\\varepsilon_x(t)$ and $\\varepsilon_y(t)$ are the pulses' amplitudes of the two independent quadrature controls (XY control)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the ideal case, we can ignore the higher energy levels of the qubit. To implement a $\\theta$ rotation about the x-axis, we set $\\delta _1$ to be zero and solve the equation directly:\n",
    "$$\n",
    "\\int_0^{t_g}\\varepsilon_x(t)dt=\\theta. \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As for a Gaussian waveform $\\varepsilon_G=Ae^{\\frac{1}{2}(\\frac{t-\\tau}{\\sigma})^2}-B$, we solve $\\int_0^{t_g}\\varepsilon_G(t)dt=\\theta_x$ to determine the amplitude $A$ corresponding to a $\\theta_x$ rotation about the x-axis:\n",
    "$$\n",
    "A=\\theta_x/\\left( \\int_0^{t_g}e^{-(t-\\tau)^2/2\\sigma^2}dt-t_ge^{-\\tau^2/2\\sigma^2} \\right),\n",
    "$$\n",
    "$$\n",
    "B=Ae^{-\\tau^2/2\\sigma^2}.\n",
    "$$\n",
    "In the equations above, $A$ ensures that the desired magnitude of rotation is implemented; while $B$ enforces that the pulse's amplitude start and end on zero. \n",
    "\n",
    "In the following code, we first define a couple of parameters to set the rotation angle and the anharmonicity of the system. Then, we define the functions for calculating parameters of the Gaussian waveform (commonly used waveform functions are available in Quanlse)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta_x = np.pi # the angle of rotation\n",
    "Delta = -0.4 * 2 * np.pi # the anharmonicity in GHz\n",
    "\n",
    "# Calculate the parameters\n",
    "def intTheta(tg):\n",
    "    y = integrate.quad(gaussian, 0, tg, {\"a\": 1, \"tau\": 0.5 * tg, \"sigma\": 0.25 * tg})\n",
    "    return y[0]\n",
    "\n",
    "def calAx(tg):\n",
    "    return theta_x / (intTheta(tg) - gaussian(0, args={\"a\": 1, \"tau\": 0.5 * tg, \"sigma\": 0.25 * tg}) * tg)\n",
    "\n",
    "def calBx(tg):\n",
    "    return calAx(tg) * gaussian(0, args={\"a\": 1, \"tau\": 0.5 * tg, \"sigma\": 0.25 * tg})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the DRAG procedure, the waveforms and detunings are modified to:\n",
    "$$\n",
    "\\varepsilon_y(t) = -\\frac{\\dot {\\varepsilon_x}(t)}{\\alpha_q}, \n",
    "$$\n",
    "$$\n",
    "\\delta_1(t) = -\\frac{\\varepsilon_x^2(t)}{2\\alpha_q}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we build the control pulses $\\varepsilon_x(t)$ and $\\varepsilon_y(t)$ and set the drive detuning $\\delta_1$ according to the equations above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the control waveforms\n",
    "def epsilonX(t, params):\n",
    "    tg = params['tg']\n",
    "    a = calAx(tg)\n",
    "    b = calBx(tg)\n",
    "    return gaussian(t, args={\"a\": a, \"tau\": 0.5 * tg, \"sigma\": 0.25 * tg}) - b\n",
    "    \n",
    "def epsilonY(t, params):\n",
    "    tg = params['tg']\n",
    "    a = calAx(tg)\n",
    "    return dragY1(t, args={\"a\": a, \"tau\": 0.5 * tg, \"sigma\": 0.25 * tg}) / Delta\n",
    "\n",
    "# Set the drive detuning  \n",
    "def delta1(t, params):\n",
    "    tg = params['tg']\n",
    "    lamda = np.sqrt(2)\n",
    "    return - epsilonX(t, {\"tg\": tg}) ** 2 / 2 / Delta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quanlse realization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quanlse stores the system's information required for simulation and optimization in the Hamiltonian dictionaries. First of all, we create an empty Hamiltonian dictionary using the function `createHam()` and add terms to it using previously defined parameters and functions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initiate the Hamiltonian dictionaries\n",
    "ham = qham.createHam(title=\"no drag\", dt=0.1, qubitNum=1, sysLevel=3)\n",
    "ham_drag = qham.createHam(title=\"drag\", dt=0.1, qubitNum=1, sysLevel=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this particular task, the system Hamiltonian can be expressed in four terms:\n",
    "$$\n",
    "\\hat H_R = \\hat H_{\\rm drift} + \\hat H_{\\rm xctrl} + \\hat H_{\\rm yctrl}+ \\hat H_{\\rm freq} ,\n",
    "$$\n",
    "where $\\hat H_{\\rm drift}= \\alpha_q\\hat a^{\\dagger}\\hat a^{\\dagger}\\hat a \\hat a/2$ represents the anharmonicity of the qubit, which is intrinsic and time-independent. We add the drift terms by calling `addDrift()`. The operator $\\hat a^{\\dagger}a^{\\dagger} \\hat a \\hat a$ is defined as `duff()` in Quanlse, which takes the system's dimension as a parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add the anharmonic terms\n",
    "qham.addDrift(ham, \"drift\", onQubits=0, matrices=duff(3), amp=Delta / 2.0)\n",
    "qham.addDrift(ham_drag, \"drift\", onQubits=0, matrices=duff(3), amp=Delta / 2.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, the control terms $\\hat H_{\\rm xctrl}=\\frac{1}{2}(\\hat a +\\hat a^{\\dagger})$, $\\hat H_{\\rm yctrl}=\\frac{i}{2}(\\hat a -\\hat a^{\\dagger})$ and $ \\hat H_{\\rm freq}=\\hat a^{\\dagger}\\hat a $ are added by calling the function `addControl()`. Their according operators are also available and can be found in `Utils.Operator`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add the control terms\n",
    "qham.addControl(ham, \"ctrlx\", onQubits=0, matrices=driveX(3))\n",
    "qham.addControl(ham_drag, \"ctrlx\", onQubits=0, matrices=driveX(3))\n",
    "qham.addControl(ham_drag, \"ctrly\", onQubits=0, matrices=driveY(3))\n",
    "\n",
    "# Add the detuning term\n",
    "qham.addControl(ham_drag, \"detune\", onQubits=0, matrices=number(3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a thorough comparison, we compute gate fidelities within a range of gate times. In fact, the task can be done very efficiently using Quanlse. In particular, `runHamiltonian()` supports batch-job simulation and returns a list of dictionaries with details of the result, and the unitary operator is stored under the key `\"unitary\"`.\n",
    "\n",
    "The simulation may take a long time to process on local devices. However, Quanlse provides a cloud service that could speed up this process significantly. To use Quanlse Cloud Service, the users can get a token from http://quantum-hub.baidu.com and use the functions in `runHamiltonian()` module to submit the job onto Quanlse's server. \n",
    "\n",
    "After the simulation, we assess the performance of the implemented gate using DRAG pulse by calculating the infidelity for various gate time defined as:\n",
    "\n",
    "$$\n",
    "{\\rm infid} =1- \\frac{1}{2}\\left|{\\rm Tr}(\\hat{\\sigma}_x P(U))\\right|.\n",
    "$$\n",
    "\n",
    "Here, the projected evolution $P(U)$ ($U$ is the evolution of the system) in particular describes the evolution projected to the computational subspace consisting of the two lowest energy eigenstates $|0\\rangle$ and $|1\\rangle$; $\\hat{\\sigma}_x$ is the target gate we want to implement. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Gate times at which to compute gate fidelities\n",
    "t = np.arange(2., 9., 0.5)\n",
    "\n",
    "# Create the arrays for storing gate infidelities\n",
    "errorx = np.zeros(len(t))\n",
    "errorxdrag = np.zeros(len(t))\n",
    "\n",
    "# Intialize array index\n",
    "jobList = []\n",
    "jobList_drag = []\n",
    "for tg in t:\n",
    "    jobWaves = []\n",
    "    jobWaves_drag = []\n",
    "    # Add Gaussian Wave of X control on the qubit 0\n",
    "    paraArgs = {\"a\": -0.5 * 2.0 * np.pi}\n",
    "    # Add wave for the job list without DRAG pulses\n",
    "    jobWaves.append(makeWaveData(ham, \"ctrlx\", f=epsilonX, para={\"tg\": tg}, t0=0, t=tg))\n",
    "    # Add wave for the job list with DRAG pulses\n",
    "    jobWaves_drag.append(makeWaveData(ham_drag, \"ctrlx\", f=epsilonX, para={\"tg\": tg}, t0=0, t=tg))\n",
    "    jobWaves_drag.append(makeWaveData(ham_drag, \"ctrly\", f=epsilonY, para={\"tg\": tg}, t0=0, t=tg))\n",
    "    jobWaves_drag.append(makeWaveData(ham_drag, \"detune\", f=delta1, para={\"tg\": tg}, t0=0, t=tg))\n",
    "    # Append this job to the job list\n",
    "    jobList.append(jobWaves)\n",
    "    jobList_drag.append(jobWaves_drag)\n",
    "\n",
    "# Submit the job lists to Quanlse Cloud Service\n",
    "result = runHamiltonian(ham, jobList=jobList)\n",
    "result_drag = runHamiltonian(ham_drag, jobList=jobList_drag)\n",
    "errorx = []\n",
    "errorx_drag = []\n",
    "for index in range(len(t)):\n",
    "    errorx.append(unitaryInfidelity(np.array([[0, 1], [1, 0]], dtype=complex), result[index][\"unitary\"], 1))\n",
    "    errorx_drag.append(unitaryInfidelity(np.array([[0, 1], [1, 0]], dtype=complex), result_drag[index][\"unitary\"], 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can analyze and visualize the results using the Matplotlib library."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.semilogy(t, errorx_drag, label='With DRAG', marker='.')\n",
    "plt.semilogy(t, errorx, label='Without DRAG', marker='.')\n",
    "\n",
    "plt.xlabel('Gate Time (ns)')\n",
    "plt.ylabel('Infidelity')\n",
    "plt.title('X Gate')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As demonstrated above, most of the leakage error is mitigated. The blue (DRAG optimized waveform) line illustrates that DRAG reduces the infidelity by orders of magnitude."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "This tutorial introduces the DRAG technique using Quanlse. The users are encouraged to explore other advanced research which are different from this tutorial."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\\[1\\] [Motzoi, Felix, et al. \"Simple pulses for elimination of leakage in weakly nonlinear qubits.\" *Physical review letters* 103.11 (2009): 110501.](https://link.aps.org/doi/10.1103/PhysRevLett.103.110501)\n",
    "\n",
    "\\[2\\] [Krantz, Philip, et al. \"A quantum engineer's guide to superconducting qubits.\" *Applied Physics Reviews* 6.2 (2019): 021318.](https://aip.scitation.org/doi/abs/10.1063/1.5089550)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
